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Determination of the local config^ration of interacting defects in a crystalline, periodic solid is problematic 
because defects typically do not have a long-range periodicity. Uranium dioxide, the primary fuel for fission 
reactors, exists in hypers toichiometric form, UO2+X. Those excess oxygen atoms occur as interstitial defects, 
and these defects are not random but rather partially ordered. The widely-accepted model to date, the Willis 
cluster based on neutron diffraction, cannot be reconciled with the first-principles molecular dynamics 
simulations present here. We demonstrate that the Willis cluster is a fair representation of the numerical ratio 
of different interstitial O atoms; however, the model does not represent the actual local configuration. The 
simulations show that the average structure of UO2+X involves a combination of defect structures including 
split di-interstitial, di-interstitial, mono-interstitial, and the Willis cluster, and the latter is a transition state 
that provides for the fast diffusion of the defect cluster. The results provide new insights in differentiating the 
average structure from the local configuration of defects in a solid and the transport properties of UO2+X. 

Uranium dioxide is the principal fuel of nuclear reactors. One of the unique properties of UO2 is its ability to 
accommodate a variable stoichiometry, depending on temperature and oxygen pressure^'^. The excess 
oxygen atoms in hyperstoichiometric uranium dioxide (U02+x) occur as interstitial defects^"^\ Positions 
and dynamics of these excess oxygen atoms control many important properties, such as thermal conductivity^^"^^, 
fission-product accommodation and transport^^, micro -structure evolution^^'^^, and corrosion behavior^^'^^. 
These properties are closely related to the performance of the fuel in a reactor and its behavior in a geologic 
disposal. For UO2+X at low x values, the interstitial O atoms occur as isolated point defects. As x increases, 
individual defects interact with each other increasingly and form clusters^"^'^°'^\ Various experimental and 
theoretical studies have shown that these clusters are not random but rather structured with well-defined con- 
figurations^"^'^ However, the defect structures are difficult to quantify experimentally using diffraction 
techniques because the clusters are local structures and lack long-range periodicity. Based on early neutron- 
diffraction studies, a defect cluster model, the so-called 2:2:2 Willis type, was proposed for UO2.1 1-2.13 over fifty 
years agol Since then, this has remained the dominant conceptual model in the literature. Although slightly 
modified later by the original author and the collaborators^'^, the proposed oxygen configuration remains the 
same'*'^. The acceptance of the model is largely based on the neutron diffraction data. However, recent first- 
principles calculations and empirical potential molecular dynamics show that the Willis cluster is not stable^ ^'^^"^^ 
Upon optimization, it spontaneously relaxes to a split di-interstitial or tri- oxygen cluster sharing a vacancy 
(V30")^^'^^"^^. The essential question is whether the Willis defect is an appropriate model for UO2+X or whether 
it is a limitation of static first-principles calculations being performed at the athermal limit that cannot account for 
finite temperature effects on the defect structure. In order to overcome the limitations related to zero temperature 
and the accuracy of empirical potentials in the previous theoretical calculations^ ^'^^"^^ first-principles molecular- 
dynamics simulations at high temperatures are employed here. These results provide a self- consistent explanation 
of the Willis cluster model that is based on neutron diffraction data and the atomistic- scale structure that is 
derived from recent theoretical calculations. The simulations also improve the understanding of different defect 
types and reveal the role of the Willis defect model in the transport of the oxygen defect clusters in U02+x- 



Results and discussion 

Uranium dioxide has the isometric fluorite structure (Fm3m). The 4a site is occupied by uranium, the 8c site by 
oxygen, and 4b site is empty in UO2. The Willis 2:2:2 defect cluster consists of two vacancies (Vo) at the O 8c site, 
two O interstitials displaced in <110> directions (O') from the 4b site, and two other O interstitials displaced in 
< 1 1 1> (O") from the 4b site (Figure la)^ Starting with this configuration without constraints, the cluster was 
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optimized to a split di-interstitial defect (V30"). Depending on how 
the WilUs defect is perturbed by a small change in geometry, one of 
the < 1 1 1 > oxygen interstitials (O") moves back to one vacant lattice 
8c site while the other O" and two O' move away from their initial 
positions towards to the other vacant lattice 8c site. This split di- 
interstitial defect (V30") consists of three interstitial O atoms shar- 
ing an oxygen 8c site. This result confirms the instability of the Willis 
defect cluster repeatedly reported in recent first-principles studies 
using similar computational methods^^"^^. Geometry optimizations 
starting from the Willis defect configuration (Figure la) and with a 
constraint to allow the two <110> O atoms to move only in the 
(001) plane result in a only slightly modified Willis cluster config- 
uration. However, when the constraint is removed, the structure 
relaxes to a split di-interstitial defect. This result suggests that, 
although the Willis cluster is unstable, it could provide a transition 
path for cluster migration from one location to another, very differ- 
ent from previously studied paths for oxygen di-interstitial dif- 
fusion^ In order to check the energy barrier of the transition 
path, the Climbing Image-Nudged Elastic Band (CI-NEB) method^^ 
was applied to estimate the barrier for the Willis cluster configuration 
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Figure 1 | Schematic structure of Willis 2:2:2 defect cluster (a), which 
relaxes to either one of the split di-interstitial defects, depending on the 
path I or II respectively. The arrows indicate the direction of the 
movements of the atoms, (b) Potential energy surface of the migration of a 
split di-interstitial defect cluster. Migration energy barrier: 0.14 eV. 
Arrows in (b) show relaxation directions. Spheres are O atoms. Uranium 
atoms are not shown. Red stands for lattice 8c O atoms, purple and pink for 
oxygen interstitials displaced along <111> <110>, respectively, and grey 
circles for the oxygen vacancy at the 8c site. 



as a transition state. The calculated migration energy is —0.14 eV 
(Figure lb), which is low, about 5 times ksT at ambient temperature. 
Similar migration energies were estimated using empirical potential 
models and temperature-accelerated molecular-dynamics simula- 
tions but assuming different transition paths^^'^^. A recent DFT cal- 
culation estimated a higher migration barrier of 0.47 eV for a split di- 
interstitial assuming a different path^^ Since static diffusion calcula- 
tions often have to presume a transition path, one may question 
whether the paths and calculated values at the athermal limit repres- 
ent the properties of the diffusion at the finite temperatures at which 
experiments with UO2+X are often performed. 

In order to reconcile the structure of the Willis defect and to 
directly reveal the migration of a split di-interstitial defect cluster, 
first-principles molecular- dynamics (MD) simulations were carried 
out at 300 K, 500 K, 800 K, 1200 K, 1600 K, and 2000 K. It was 
expected that thermal energy at higher temperatures activates the 
cluster migration and increases the effectiveness of probing the phase 
space of the cluster migration, giving its low migration energy bar- 
rier. The simulations started with one split di-interstitial in the com- 
putational supercell (UO2.06)- At low temperatures (300 K-800 K), 
the split di-interstitials remain at their initial location with a config- 
uration of the two interstitials displaced from the octahedral 4b sites, 
pushing a lattice 8c oxygen atom towards the other octahedral 4b site, 
resulting in three interstitials sharing an oxygen lattice 8c site. 
Figure 2a shows the atomic density probability contour maps at 
300 K projected on (001) plane. Similar results were observed at 
500 K and 800 K. At 1200 K, one of the interstitial O atoms is moved 
back to its lattice site while another oxygen interstitial appears, 
resulting in a new split di-interstitial defect cluster at a different 
location (Figure 2b). At 1600 K, the cluster has become two separate 
mono -interstitials (Figure 2c). After the dissociation, the two mono- 
interstitial defects remain relatively immobile for the rest of the 
molecular dynamics simulation. As the temperature increases to 
2000 K, the cluster becomes highly mobile and goes through mul- 
tiple transformations at different locations (Figure 2d). 

A careful analysis of the trajectory of UO2.06 at 2000 K shows that, 
as the split di-interstitial cluster migrates, it is transformed from one 
to another of the same kind. During the transition, the defect cluster 
passes quickly through a transition state with two 8c site vacancies 
(Vo), two <110> interstitials (O'), and two <111> interstitials 
(O"), which resembles the 2:2:2 Willis defect cluster. Between the 
transitions, the split di-interstitials cluster remain relatively stable. At 
some point during the MD run, the split di-interstitial cluster was 
temporarily transformed into a di-interstitial (i.e., two mono-inter- 
stitials at neighboring locations), but transformed back to a split 
di-interstitial cluster, which does not involve the 2:2:2 Willis con- 
figuration. Figure 3 shows the MD profiles of the potential energy, 
temperature, and pressure with time as the cluster migrates. Split di- 
interstitial clusters constitute the majority of the simulation time, and 
the Willis cluster is only a transition state occurring in a small 
percentage of the time. In order to check how individual split di- 
interstitials clusters interact with each other, similar molecular- 
dynamics simulations of UO2.13 and UO2.19 were completed at 
2000 K with two and three di-interstitials in a 2X2X2 supercell, 
respectively. At the beginning of the simulations, each individual 
cluster behaved similarly to a single split di-interstitial cluster. 
However, once they were in contact and aggregated, the mobility 
of the clusters slowed because of complex interactions between them. 
During the migration of the interacting clusters, one of the split di- 
interstitial clusters had to be activated to a transition state similar to 
the Willis cluster. One important result of the simulations of single 
and multiple di-interstitial clusters is that the 2:2:2 Willis defect 
configuration is actually a transition state for migrating clusters at all 
of the X values investigated. 

While the oxygen defect clusters migrate, the uranium atoms are 
immobile. With the excess oxygen atoms in the structure, some 
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Figure 2 | Atomic density probability contour maps of UO2.06 (U32O66) projected onto (100) plane, calculated from the MD trajectories at different 
temperatures (a: 300 K, b: 1200 K, c: 1600 K, d: 2000 K). Circles indicate the locations of the interstitial O atoms. 



uranium atoms are oxidized from U^^ to U^^ by transferring an 
electron from uranium to oxygen. An interesting observation is that, 
at low temperatures (300 K-800 K), the U^^ atoms are mostly adja- 
cent to the defect cluster, while at temperatures above 1200 K, the 5 + 
oxidation state occurs at all uranium lattice sites over the course of 
the simulation time of the 3 ps as shown in Figure SI . The fact that the 
location of U^^ is not constrained by the interstitial O atoms at 
temperatures lower than the temperatures at which the di-interstitial 
cluster becomes mobile suggests that the charge mobility is 
decoupled from the interstitial oxygen mobility. It seems that the 
oxidation state of U atoms changes between U^^ and U^^ at a faster 
rate than oxygen cluster migration through the U lattice. 

In order to reconcile these computational results with the early 
neutron diffraction data, average occupancy numbers and displace- 
ments of O' and O" interstitials were calculated over the trajectory for 
each composition and are listed in Table 1. For UO2.06, at low tem- 
peratures, the interstitials are mainly displaced along < 111 >. As the 
temperature increases, the 070" ratio increases. As expected, the 
displacements increase with temperature as well. As x increases from 
2.06, 2.13, to 2.19 at 2000 K, the 070" ratio increases because the 
interstitial O atoms joining neighboring di-interstitial clusters are 
mainly displaced along <110>. The calculated ratios of O' and O" 
occupancies for UO2.13, 0.15:0.09, is consistent with the experi- 
mental values of 0.08-0.33:0.10-0.16 for UO2.11-UO2.13 based on 
neutron diffraction^ The calculated displacements for UO2.13 are 
0.77 ± 0.21 A and 0.92 ± 0.20 A along <110> and <111>, respect- 
ively, as compared with experimental values of 0.85 ± 0.08 A and 
1.04 ± 0.10 A for UO2.12 at 1073 Note that the calculations were 
done at a higher temperature for the composition, which is necessary 
to have adequate statistical averages in a short MD simulation. 
However, the temperature only has a small effect on the values of 
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Figure 3 | Potential energy, temperature, and pressure profiles over the 
course of a —3.5 ps MD simulation of UO2.06 (U32O66) at 2000 K. Top 

panel shows defect cluster configurations of different types: dsi split di- 
interstitial (blue), w: Willis 2:2:2 defect cluster (red), d: di-interstitial 
(green). 
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Table 1 Average occupancy numbers and displacements of <1 1 0> and <1 1 1 > interstitial O atoms in 
different temperatures, v is the displacement in < 1 1 0> directions and w is in the < 1 1 1 > directions 


UO2+X with different x values at 


Composition 




UO2.06 (U32O66) 




UO2.13 


UO2.19 


Temperature 


300 K 


1200 K 


2000 K 


2000 K 


2000 K 


O/O'/O" 
v<110> 
v/<lll>(A) 


1 .97/0.00/0.09 
0.67 ± 0.07 
0.78 ± 0.09 


1 .95/0.05/0.06 
0.71 ± 0.18 
0.87 ±0.19 


1 .95/0.06/0.05 
0.74 ± 0.27 
0.90 ± 0.25 


1.89/0.15/0.09 
0.77 ±0.21 
0.92 ± 0.20 


1.90/0.19/0.10 
0.72 ± 0.24 
0.90 ± 0.24 



both the 070" ratio and displacements as the temperature increases 
from 1200 K to 2000 K as shown in Table 1 for UO2.06. The lower 
temperature is comparable to the experimental temperature con- 
ducted at 1073 K^. The probability distributions of the angle between 
< 1 1 1 > and the displacement direction of interstitial oxygen atom 
are shown in Figure S2. The calculated result suggests that the 2:2:2 
Willis cluster model for U02.ii_i3 does account for the numerical 
fraction of the O' and O" inter stitials. However, the Willis model 
does not represent the local defect configuration of defect clusters. 
The often assumed Willis 2:2:2 defect configuration is, in fact, a 
transition state for the migration of a split di- interstitial cluster in the 
hypers toichiometric UO2. A careful review of all the trajectories 
shows that the average structure of UO2+X involves a combination 
of defect structures including the split di-interstitial, di-interstitial, 
mono -interstitial, and Willis cluster, and the latter serves as a trans- 
ition state for a fast diffusion of the defect cluster (Figure 4). 
Spectroscopic techniques such as vibrational spectroscopy can be 
used to validate local structure configuration of different types of 
defect clusters, and property measurements such as electron and 
ionic conductivity can be used to test if the charge transport between 
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Figure 4 | Defect models for oxygen interstitials in UO 2+x- The Willis 
defect cluster (a) serves as a transition state for rapid diffusion of the split 
di-interstitial defect (b). The latter can also migrate through a di- 
interstitial (c). Di-interstitial can dissociate to two immobile mono- 
interstitials (d). 



U^^ and U^^ is activated at a lower temperature than the oxygen 
migration. 

Methods 

The first-principles calculations were based on the Density Functional Theory and 
plane wave basis sets as implemented in VASP^^. The Projector-Augmented Wave 
method and exchange-correlation as parameterized by the Perdew-Wang 91 func- 
tional were applied in the Generalized Gradient Approximation^*^-^''. For U, 14 elec- 
trons are treated as valence electrons and the core electrons have [Xe, 5d, 4f] 
configuration. For O, 6 electrons are treated as valence electrons and the core has 
[He]. On-site Coulomb interaction (U=4.5, J=0.51) with the rotational invariant 
Liechtenstein approach^", fully relativistic calculation for the core-electrons, and the 
scalar relativistic approximation for the valence electrons were employed to account 
for electron correlation and relativistic effects^\ The U and J values are based on the 
experiments^, and these values have been used in a number of most recent publica- 
tions of DFT-HU studies of 1102*^'^^'^"^. However, other different U and J values have 
also been used in the literature^^'^^, depending on how U and J values were deter- 
mined. We have used U=3.8 eV and J=0.4 eV in previous publicationsS'^'S^, the 
results are comparable with those using U = 4.5 eV and J = 0.51 eV. In this study, for 
the reason of comparison and consistency with the majority of the literature, the 
values based on the experiment were used. The details and performance of the 
methods on the calculations of crystal structure and electronic structure of UO2+X are 
summarized in a previous publication^^. All the calculations were performed with a 
supercell of 2X2X2 unit cells of UO2 and the energy cut-off for the plane-wave basis 
was set to 520.00 eV. 3X3X3 k-point grids were used for the static calculations and T 
point for the molecular dynamics simulations. Ferromagnetic ordering without spin- 
orbit coupling was used in the molecular-dynamics simulations. Antiferromagnetic 
and ferromagnetic configurations were both used in the static calculations to test the 
effect of the spin configurations on the calculations. The optimized volumes were 
used for static calculations. For dynamics simulations, both optimized volumes and 
relaxed volumes at 1 bar and high temperatures were used. Molecular-dynamics 
simulations were performed using a NVT ensemble with the Nose-thermostat for 2-3 
ps for equilibration, which was checked by monitoring a number of parameters 
including temperature, pressure, potential energy of the systems, and kinetic energy 
of both U and O sublattices. The time step is 0.5 fs for low- temperature simulations 
(300-800 K) and 0.25 fs for high-temperature (1200-2000 K). It is necessary to carry 
out the simulations at high temperatures, at which multiple migration events can be 
observed and statistical meaningful averages can be obtained in short molecular 
dynamics simulations. The total energy drift was ~ (2-8) meV/atom/ps at low 
temperatures and increased to —25 (15-50) meV/atom/ps at high temperatures. For 
UO2+X with unpaired 5f electrons, it was extremely challenging to achieve both high 
accuracy and adequate sampling with enough statistics of the systems. Extensive tests 
were performed to balance between the accuracy and duration of the simulations by 
systematically tuning simulation parameters including those related to convergences, 
real space projection, and time step. After equilibrations, simulations run up to 4 ps 
for statistical analysis. 
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